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Two-dimensional crystals are emerging materials for nanoelectronics. Development of the 
field requires candidate systems with both a high carrier mobility and, in contrast to 
graphene, a sufficiently large electronic bandgap. Here we present a detailed theoretical 
investigation of the atomic and electronic structure of few-layer black phosphorus (BP) to 
predict its electrical and optical properties. This system has a direct bandgap, tunable from 
1.51 eV for a monolayer to 0.59 eV for a five-layer sample. We predict that the mobilities are 
hole-dominated, rather high and highly anisotropic. The monolayer is exceptional in having 
an extremely high hole mobility (of order 10,000 cm^V~^ s~b and anomalous elastic 
properties which reverse the anisotropy. Light absorption spectra indicate linear dichroism 
between perpendicular in-plane directions, which allows optical determination of the 
crystalline orientation and optical activation of the anisotropic transport properties. These 
results make few-layer BP a promising candidate for future electronics. 
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The discovery of graphene laid the foundations for 
many new areas of research. One of the most 
important foundation is the intensive investigation of 
two-dimensional (2D) atomic-layer systems, including graphene 
itself^"^, transition metal dichalcogenides (TMDCs)^"^, 
silicene^'^^ and germanane^^'^^ as candidate materials for 
future electronics appHcations^^"^^. A high-performance device 
such as a field-effect transistor (FET) requires a moderate 
electronic band gap, a reasonably high carrier mobility of 
the channel material and excellent electrode- channel 
contacts"*'^'^'^^"^^. Graphene offers extremely high mobilities, 
due to its very low carrier effective mass, and thus is considered to 
be a promising candidate for high-speed FET devices, but its 
intrinsic dispersion is gapless^""*' . Despite extensive efforts 
following a wide variety of approaches to the problem of opening 
a gap in different graphene nanostructures, all devices to date 
have a relatively large 'off current and thus a low 'on-off ratio^^. 
For this reason the emergence of monolayer TMDs has attracted 
substantial research interest, and in particular M0S2 has recently 
been used to fabricate a FET^^. Unlike graphene, mono- 
layer M0S2 is a direct-bandgap semiconductor with a carrier 
mobility of ~200 cm^ V~ ^ s~ ^ (ref. 15), improvable up to 
500cm V~ ^ s~ ^ (ref. 8), a value which is not unreasonable for 
applications but remains orders of magnitude lower than that of 
graphene^^'^"*. Thus the entire community working on atomic- 
layer transport continues to search for a 2D material which is 
semiconducting, preferably with a direct gap, has high carrier 
mobility and has the potential to form excellent contacts with 
known electrode materials. 

Here we present the theoretical discovery of a new category of 
layered direct-bandgap semiconductor, few-layer black phos- 
phorus (BP), with high mobility, high in-plane anisotropy and 
linear dichroism. BP is an allotrope of phosphorus, whose crystal 
structure is a strongly folded honeycomb sheet with 'troughs' 
running along the y-axis (b direction). We demonstrate using 
density functional theory (DFT) calculations that few-layer BP 
systems, from the monolayer up to five-layer structures, are 
thermally stable with an interlayer interaction energy of 
— 0.46 eV and a bandgap that falls exponentially with the 
thickness. Carrier mobilities at room temperature are high, and 
furthermore exhibit strongly anisotropic behaviour. The BP 
monolayer has an exceptionally high hole mobility, which we 
predict to lie in the range from 10,000 to 26,000 cm^ V ~ ^ s ~ ^ 
and which we explain by visualizing the wavefunctions. We find 
an expHcit linear dichroism in the computed absorption spectra, 
whereby the positions of the lowest-energy absorption peaks 
for the two in-plane directions differ strongly. These results 



demonstrate that few-layer BP is a new category of 2D 
semiconductor with high potential for novel applications in 
nanoelectronics and optoelectronics. 

Results 

Geometric and electronic properties of bulk BP. Because no 
experimental data are available for the atomic structure and 
electronic properties of few-layer BP, we begin by considering 
bulk BP to refine the accuracy of our theoretical predictions. We 
have performed DFT calculations using a number of different 
functional to gauge which provide the best fit to experiment; a 
full list may be found in Supplementary Table 1 and we sum- 
marize our findings here. If one considers only the lattice geo- 
metry then the PBE-G06 (refs 20,21) and optB86b-vdW^^ 
methods produce the best results. However, neither local 
density approximation-modified Becke-Johnson (LDA- 
j^gj^23,24 j^Qj. HSE06 (refs 25,26) calculations based on these 
'best' geometries can provide a satisfactory value for the bulk 
bandgap, as shown in Supplementary Fig. 1. To optimize both 
sets of properties simultaneously, we found that the optB88-vdW 
functional combined with either the LDA-mBJ or HSE06 
method to predict the electronic bandstructure, gave the best fits. 
Details of choosing the best-fit functional are available in 
Supplementary Methods. Figure la shows the fully relaxed 
atomic structure of bulk BP, where the equilibrium lattice 
constants obtained from optB88-vdW are only 1-2% larger than 
the experiment^^"^^. The accompanying Brillouin zone (BZ) and 
electronic bandstructures, labelled mBJ (optB88-vdW) and 
HSE06 (optB88-vdW) are shown in Fig. lb,c. Both 
combinations of methods predict that bulk BP is a 
semiconductor with a direct bandgap at the Z point of 0.31 eV 
(mBJ) or 0.36 eV (HSE06), both values being fully consistent with 
the experimental value of 0.31-0.35 eV (refs 31-34). 

At the bandgap wavevector, which is the Z point, one valence 
band (VB) and one conduction band (CB) disperse rather 
strongly along the Z-Q and Z-G directions (Fig. Ic), indicating 
very small effective masses. A fit of these bands using the nearly- 
free electron model gives effective carrier masses along the Z-Q 
direction which are rather small and similar, namely 0.12 mo for 
electrons and 0.11 mo for holes, whereas these along Z-G are 
slightly larger, taking the respective values 0.15 mo and 0.30 mo 
(all masses taken from HSE06). Considerably larger values are 
found along Z-T'-A^ where the carrier effective masses are 1.15 
mo and 0.71 mo, respectively. All of these results are fully 
consistent with the experimental values for bulk BP^^. The 
corresponding mBJ results are within 0.04 mo of the HSE06 
values along Z-Q, Z-G and Z-T'-A' except for the electron 
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Figure 1 | Lattice and electronic structures of bulk black phosphorus, (a) Crystal structure of bulk BP marked with coordinate axes (x, y, z), lattice vectors 
(0, b, c) and structural parameters {R\, R2, ^^ and 92). (b) Brillouin zone path of BP primitive cell, (c) Electronic bandstructures for bulk BP calculated 
with the HSE06 functional (red solid line) and the mBJ potential (blue dashed line), together with fitted effective masses along the Z-T'-A', Z-Q 
and Z-G directions. At the right of the image, a zoomed-in plot shows the direct bandgap at Z. ^vbm is the energy of valence-band maximum. 
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effective mass along Z-T'-A^ which is larger than the HSE06 
value by 0.13 mo- Values for the direction Z-Q are similar to 
other high-mobility semiconductors, such as AlGaAs (0.1 mo)^^, 
and four times smaller than those of M0S2 (ref. 37). 

Geometric and electronic properties of few-layer BP. The 

similarity between the bandstructures of bulk and monolayer 
BP^^, and the low effective masses of bulk BP^^, suggest that few- 
layer BP is likely to be a high -mobility, direct-bandgap 2D 
semiconductor. We have thus performed direct calculations of the 
geometric and electronic properties for few-layer BP systems. 
Table 1 summarizes the changes of geometrical properties as a 
function of the layer number from 1 to 5. The lattice parameter a 
increases by 0.11 A on passing from bulk to monolayer BP, 
whereas b grows only by 0.02 A. There is an abrupt reduction of a 
between the monolayer and the bilayer, which we believe to be the 
primary consequence of interlayer interactions (0.46 eV) in the 
bilayer. The significant stretching of a in few-layer BP systems is 
caused almost exclusively by an increase of the bond angle 61 but 
not of the bond lengths. 

Results from our bandstructure calculations for the five few- 
layer BP systems are shown in Fig. 2. In the monolayer, the 
original Z point of the bulk BZ folds back to the Gamma (G) 
point, so that the original Z-Q and Z-T'-A' directions of the bulk 
BZ project onto the G-X and G-Y directions of the monolayer, 
which correspond respectively to the a and b directions in real 
space (Fig. 2a,b). Results obtained from the mBJ method are 
quantitatively the same as those obtained from HSE06. 
Supplementary Figure 2 and Supplementary Table 2 show that 
the predicted bandgap is rather less sensitive than the unit- cell 
dimensions to the functional used for optimizing the atomic 
structures, with similar values emerging for all few-layer BP 
systems from the best-performing functional. However, other 
physical quantities, notably the location of band maximum 
(minimum) and the carrier effective mass (and hence the 
mobility) are in fact very sensitive to the choice of the functional, 
and hence we have focused primarily on our optB88-vdW results 
for the discussion of these. Monolayer BP is indeed a direct- 
bandgap semiconductor, as shown in Fig. 2c, with the gap value 
of 1.51 eV obtained at the G point. Phosphorus is a heavier 
element than carbon and therefore should have stronger spin- 
orbit coupling (SOC) in its 2D forms than graphene. We have 
considered full SOC effects in calculating the electronic 
bandstructure of bulk BP. Supplementary Figure 3 shows that 
inclusion of SOC terms has no appreciable effect on the primary 
features of the bandstructure, indicating that phosphorus is 'not 
sufficiently heavy' to cause any quaHtative changes. Only a very 
small separation of the formerly fourfold degenerate bands, 
into two band pairs, can be found from X along N (or M for 
few-layer BP) to Y, which reaches a maximum of 23 meV around 
the N (M) point. Along G-X and G-Y, these bands are already 



separated by bandstructure effects and SOC causes no further 
splitting. That the spin- orbit splitting energy is so small is a 
consequence of the small effective nuclear charge (Zgff) of the P 
atom and the weak variation of the charge gradient in an 
elemental system such as BP. 

When two monolayers are combined to form a bilayer 
(Fig. 2d,e), the gap is reduced to 1.02 eV and two additional 
bands emerge around the gap at the G point. Together with the 
original VB and CB, we denote these bands as VBl, VB2, CBl 
and CB2, as shown in Fig. 2f. In real space, the states from these 
four bands with k close to G are extended throughout the bilayer, 
as shown in Fig. 2g. States VBl and VB2 differ in the interlayer 
region although they share the same origin in terms of the atomic 
orbitals. A clear bonding-like feature is visible in the interlayer 
region (marked by red rectangles) for VBl, which lies lower in 
energy, whereas VB2 shows an anti-bonding feature. Similar 
bonding and anti-bonding features are also found in CBl and 
CB2, where they are observed not between the layers but across 
the troughs. These features indicate that wavefunction overlap, 
rather than van der Waals effects, plays the primary role in 
mediating the interlayer interaction, which explains the abrupt 
reduction of the lattice constant a from monolayer to bilayer. 
This interlayer interaction introduces band dispersions of VBs 
and CBs in the direction normal to the layers, which leads to the 
reduction of the bandgap (by 0.5 eV) from monolayer to bilayer. 

The thicker the few-layer BP system, the stronger is the 
interlayer interaction and thus the larger is the dispersion of VBs 
(CBs), resulting in a smaller bandgap. Thus the gap falls 
continuously on adding more layers to the bilayer, reaching 
0.59 eV in five-layer BP, as illustrated in Fig. 2h. We fitted these 
values by an exponential decay relation and found that the 
corresponding bulk gap obtained by extrapolation is 0.53 eV, a 
value 0.17 eV larger than the one we calculated for bulk BP. We 
ascribe this difference to the elongated lattice parameter a in few- 
layer BP and results calculated with a series of constrained a 
values confirm our expectation; these results imply that the 
bandgap is very sensitive to the lateral strain along the a (x) 
direction and could be modulated by varying interatomic 
separations and angles. During the preparation of this manu- 
script, we became aware that a strain- induced bandgap modula- 
tion has recently been demonstrated theoretically^^. 

The small effective masses remain in all of the few-layer BP 
systems (Table 2). For the G-X direction in the monolayer, 
equivalent to Z-Q in the bulk, the carrier effective masses are 0.15 
mo (hole) and 0.17 mo (electron), only 0.04-0.05 mo larger than 
those in bulk BP. It is remarkable that in the G-Y direction, 
equivalent to Z-T'-A' in the bulk, the VB appears to be nearly flat 
close to the G point, with an effective mass of 6.35 mo (see 
Fig. 2c), nine times its bulk value of 0.71 mo. By contrast, the 
effective mass for the CB is 1.12 mo, very close to its bulk value of 
1.15 mo. Unlike the bandgap, only the hole effective mass along 



Table 1 | Structural Information for few-layer BP. 



Nl 


a (A) 


b(A) 


Ac (A) 


/?l (A) 


R2 (A) 


ever n 


e2/e2' n 


1 


4.58 


3.32 


3.20 


2.28 


2.24 


103.51 


96.00 


2 


4.52 


3.33 


3.20 


2.28 


2.24 


102.96 


96.21/95.92 


3 


4.51 


3.33 


3.20 


2.28 


2.24 


102.81/102.74 


96.30/95.99 


4 


4.50 


3.34 


3.20/3.21 


2.28 


2.24 


102.76/102.67 


96.34/96.01 


5 


4.49 


3.34 


3.20/3.21 


2.28 


2.24 


102.71/102.63 


96.37/96.05 


Bulk 


4.47 


3.34 


3.20 


2.28 


2.25 


102.42 


96.16 



BP, black phosphorus. 

Lattice constants o, b and Ac (the interlayer spacing between two adjacent BP layers) and structural parameters Rl, R2, 91 and 62 of few-layer and bulk BP calculated using optB88-vdW. There are slight 
differences in 61 and 62 between the outermost and inner layers. 
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Figure 2 | Electronic structures of few-layer BP. (a,b) Top view of the atomic structure of the monolayer and the associated Brillouin zone. (d,e) Side views 
of the atomic structure of the bilayer. (c,f) Bandstructures of monolayer and bilayer BP calculated with the HSE06 functional (red solid lines) and the mBJ 
potential (blue dashed lines), respectively. Two valence (VB1 and VB2) and two conduction states (CB1 and CB2) are marked in panel f. (g) Spatial 
structure of wavefunctions for the four marked states illustrated in the xz and yz planes using an isosurface of 0.0025 e A (h) Evolution of the direct 
bandgaps as a function of the sample thickness. Functionals used for structural optimization are shown in parentheses. Bandgap values are marked 
for the monolayer system, for the extrapolation of our results and for real bulk BP. 



Table 2 | Predicted carrier mobility. 



Carrier type 


Nl 


niyJmQ 




fix 




2D 


Cy 2D 


Hx_2D 








G-X 


G-Y 


(eV) 




(Jir 




(lO^cm^V 


-^s-i) 


e 


1 


0.17 


1.12 


2.72 ±0.02 


7.11 ±0.02 


28.94 


101.60 


1.10-1.14 


-0.08 




2 


0.18 


1.13 


5.02 ±0.02 


7.35 ±0.16 


57.48 


194.62 


-0.60 


0.14-0.16 




3 


0.16 


1.15 


5.85 ±0.09 


7.63 ±0.18 


85.86 


287.20 


0.76-0.80 


0.20-0.22 




4 


0.16 


1.16 


5.92 ±0.18 


7.58 ±0.13 


114.66 


379.58 


0.96-1.08 


0.26-0.30 




5 


0.15 


1.18 


5.79 ±0.22 


7.35 ±0.26 


146.58 


479.82 


1.36-1.58 


0.36-0.40 


h 


1 


0.15 


6.35 


2.50 ±0.06 


0.15 ±0.03 


28.94 


101.60 


0.64-0.70 


10-26 




2 


0.15 


1.81 


2.45 ±0.05 


1.63 ±0.16 


57.48 


194.62 


2.6-2.8 


1.3-2.2 




3 


0.15 


1.12 


2.49 ±0.12 


2.24 ±0.18 


85.86 


287.20 


4.4-5.2 


2.2-3.2 




4 


0.14 


0.97 


3.16 ±0.12 


2.79 ±0.13 


114.66 


379.58 


4.4-5.2 


2.6-3.2 




5 


0.14 


0.89 


3.40 ±0.25 


2.97 ±0.18 


146.58 


479.82 


4.8-6.4 


3.0-4.6 



2D, two-dimensional. 

Carrier types 'e' and 'h' denote 'electron' and 'hole', respectively. Nl represents the number of layers, nix* and my* are carrier effective masses for directions x and y, respectively, Eix (Ey) and Cx 2d (Cy_2D) 
are the deformation potential and 2D elastic modulus for the x (y) direction. Mobilities ii^_2d and iiy_2D were calculated using equation (1) with the temperature 7set to 300 K. 
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G-Y shows a strongly layer- dependent evolution, in that it 
decreases from 6.35 mo for the monolayer to 1.81 mo for the 
bilayer and eventually to 0.89 mo in the five-layer system, close to 
the bulk value of 0.71 mo. 

Carrier mobility. The electronic properties of few-layer BP are 
governed to a large extent by the carrier mobilities, which are in 
turn strongly influenced (but not solely determined) by their 
effective masses. Thus we also provide theoretical predictions for 
the carrier mobilities, in both x and y directions, in few-layer BP 
systems by applying a standard 2D model as discussed under 
Methods below. We apply a phonon-limited scattering model, in 
which the primary mechanism limiting carrier mobility is scat- 
tering due to phonons^^'"*^""*^. Two properties of the few-layer BP 
lattice, namely the deformation potential Ei and the elastic 
modulus C2D in the propagation direction of the longitudinal 
acoustic wave, are then the most relevant factors determining the 
mobility^^'^^. These quantities were calculated using the optB88- 
vdW functional and the data are presented in Table 2. If few-layer 
BP is cut into ribbons, a one-dimensional model should be 
adopted for predicting the carrier mobilities. We present the 
results and discussion for a 10-nm-width BP ribbon in the 
Supplementary Discussion and Supplementary Table 3. 

Our predicted mobilities for the five few-layer BP systems 
(Table 2) are in general moderately large (hundreds to thousands 
of cm^V~^s~^), moderately anisotropic and asymmetric 
between electrons and holes, with the holes being more mobile 
in both directions. Concerning directional anisotropy, electron 
and hole mobilities for the x direction are nearly twice the values 
along y for holes and four times as large for electrons. However, 
the monolayer is exceptional in every way. The electron mobility 
along X is nearly 14, rather than four, times the value along;/, that 
is, 1,100-1,140 cm^V~^s~^ compared with -80cm^V~^s~^ 
which is consistent with the trend of the other few-layer systems. 
To our great surprise, however, the hole mobility along x is 16-38 
times smaller than that along 7, that is, 640-700 cm^ V ~ ^ s ~ ^ 
versus 10,000-26,000 cm^V~ ^ s~ ^ making y the direction of 
higher hole conductivity. This extraordinarily large value for 
monolayer BP is a consequence of the extremely small 
deformation potential, £iy = 0.15 ± 0.03 eV, and occurs despite 
the fact that the carriers are very heavy (6.35 mo). This 
'monolayer exception' should be valuable in applications such 
as the separation of electrons and holes. 



The value of £iy for holes in a monolayer is very striking in that 
it is an order of magnitude smaller than typical values of £1, 
which are 5.0 eV for graphene"*"*, 3.9 eV for M0S2 (ref. 37), and 
3.7 eV for h-BN^^ and 7 eV for AlGaAs quantum wells^^. In few- 
layer BP systems, Eiy for holes has a conventional value 
2.97 ± 0.18 eV in a five-layer structure, decreasing smoothly to 
1.63 ± 0.16eV in a bilayer. These results can be explained rather 
well by the form of the VB wavefunctions shown in Fig. 2g, where 
interlayer and stacking- induced intra-layer overlap increases with 
layer number. These effects are completely absent in monolayer 
BP, making the VB wavefunction quite isolated along y, as it is 
also for the VB2 (anti-bonding) state of a bilayer highlighted by 
the green rectangle in Fig. 2g. Thus small structural deformations 
due to longitudinal phonon oscillations along y have almost no 
effect on this wavefunction and cause little change in its energy, 
resulting in the very small deformation potential. The situation is 
quite different for the CB (electron) wavefijnction, which has 
= 2.72 ± 0.02 eV in the monolayer and 5.02 ± 0.02 eV 
in the bilayer. The in-plane bonding and anti-bonding features 
of the CB wavefunctions in few-layer BP systems cause a 
substantial enhancement of the energy-level perturbations due to 
structural deformations along y. The monolayer exception aside, 
all El values for CB states are rather large (5.0-7.7 eV), for both x 
and y directions, and are almost independent of the number of 
layers, a result which can be explained by the strong CB 
wavefunction overlap, especially along y (Fig. 2g). For the VB 
states, the potentials are somewhat smaller (2-3 eV), 
precisely because their wavefunctions have less overlap between 
different P atoms. 

Optical absorption spectra and linear dichroism. We have also 
predicted the optical absorption spectra of few-layer BP systems 
by computing the dielectric function. Two absorption spectra are 
shown in Fig. 3a,b for light incident along the z direction and 
linearly polarized in the x and y directions, respectively; while the 
results for incident light polarized in the z direction may be found 
in Supplementary Fig. 4. These results demonstrate a strong 
linear dichroism: for a dielectric polarization in the x direction, 
the band edge of the first absorption peak is found at the bandgap 
and thus falls rapidly with the thickness of the sample (Fig. 3a). 
By contrast, with ;/-polarized light this peak is found at 3.14 eV in 
the monolayer and its position falls only slightly with thickness, 
remaining at 2.76 eV in bulk BP (Fig. 3b). From the symmetries of 



a 



< 




0.8 1.2 1.6 
Energy of incident light (eV) 



2.0 
0.0 



5.0 
0.0 



10.0 
0.0 



1 -layer 



2-layer 



5-layer 



-Bulk 



3.14 




2.90 



2.76 



5.0 
0.0 



5.0 
0.0 



2.0 



0.0 



2.0 2.4 2.8 3.2 3.6 
Energy of incident light (eV) 




Figure 3 | Optical absorption spectra. (a,b) Optical absorption spectra of few-layer BP for light incident in the c (z) direction and polarized along the 0 (x) 
and b (y) directions, respectively. Black dashed lines show an approximate linear fit used to estinnate the band edges for the first absorption peak, which are 
highly anisotropic between x and y. The band edge drops rapidly with sample thickness for x (from 1.55 eV for the monolayer (red) to 0.60 eV for five 
layers (orange)) but only slightly for y (from 3.14 eV to 2.90 eV). (c) Schematic illustration of a proposed experimental geometry to determine the orientation 
of few-layer BP structures using optical absorption spectroscopy, and thus to utilize the anisotropic electronic properties of BP. The light is linearly 
polarized in a chosen orientation and near-normally incident on the sample. The sample should be rotated to identify the a or b direction by monitoring the 
absorption signal, after which the source and drain electrodes, denoted as 'S' and 'D' respectively, may be deposited to fabricate an FET device. 
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the wavefunctions in Fig. 2g, clearly the (odd) dipole operator 
connects the VB and CB states for x-polarization, allowing the 
direct-bandgap process, but this is symmetry-forbidden for y- 
polarization and the transition occurs between VB and CB states 
elsewhere in the BZ. 



Discussion 

High carrier mobiHties would be expected in BP systems due to 
the small effective masses of 0.1-0.2 niQ. However, our 
calculations reveal that the key issue is not the effective mass 
but the very small elastic modulus C2D: all the BP systems we 
have studied are Very soft' in the x direction, due to the fact that 
bond angles are changed easily (with Httle effect on bond lengths). 
When compared to other monolayer systems, the modulus of BP 
for the X direction is 29.0 J m ~ ^, a value respectively 2, 10 and 1 1 
times smaller than M0S2 (ref 37), h-BN"*^ and graphene^^. These 
results indicate that much higher mobilities could be most likely 
achieved if few-layer BP systems could be modified to become 
stiffer, that is, to have a larger elastic modulus for the x direction, 
for example by external stress, doping or substrate confinement. 
However, possible changes of the deformation potential as a 
result of modif)^ing the elastic moduli are not yet clear and require 
further investigation. The modulus for the y direction is 
significantly larger than that for x, but the higher carrier 
effective masses, result in carrier mobilities being smaller than 
those along x. Both the bandgap and the mobility depend strongly 
on the structural properties and this offers a number of 
possibilities for further control over the electrical and 
mechanical properties of few-layer BP. 

Further, the linear dichroism we find allows an electronic 
determination of the sample orientation. We propose an 
experimental setup to deduce the orientation of few-layer 
BP systems using optical spectroscopy, as illustrated in Fig. 3c. 
If the incident light is linearly polarized in a chosen orientation 
and is near-normally incident, the absorption spectrum 
should vary when the sample is rotated and such that 
the a and b directions can be identified by monitoring the 
absorption signal. Once the sample orientation is determined, 
it is far easier to fabricate electrodes or gates utilizing the 
highest-mobility direction of a few-layer BP sample in an 
FET-type device. 

In summary, we have shown theoretically that few-layer BP is a 
novel category of 2D semiconductor offering a direct bandgap, 
high carrier mobility and high transport anisotropy, all of which 
are tunable by controlling the layer thickness. The bandgap 
decreases from 1.51 eV for a monolayer to 0.59 eV for five-layer 
BP, a set of values fitting well into the gap between graphene 
nanoribbons and TMDCs. The anisotropy of electric and optical 
properties is another unique feature distinguishing BP from other 
2D materials, such as silicene^'^^, h-BN^'^ and germanane^^'^^ 
Higher conductivity is generally found in the direction oriented 
perpendicular to the troughs (the x direction) and holes are 
usually more mobile than electrons. For holes, the mobility in the 
X direction increases from 600 cm^ V ~ ^ s ~ ^ for a monolayer to 
over 4,000 cm^ V ~ ^ s ~ ^ for five-layer BP. We found that 
monolayer BP is a very special material, in which we predict an 
extremely high hole mobility of 10,000-26,000 cm^ V ~ ^ s ~ ^ All 
of our mobility results can be understood rather well from the 
form of the real-space wavefunctions, which show a much 
stronger interlayer coupling in few-layer BP than in graphene or 
TMDCs. All these results make few-layer BP a very promising 
candidate material for future applications in electronics and 
optoelectronics. 

Note added in proof: during the review process for this 
manuscript, we became aware that few-layer BP samples have 



very recently been exfoliated from bulk BP and fabricated 
into electronic devices"*^'"*^. Measurements on these samples 
verify the high hole mobilities'*^''*^ and transport anisotropy"*^ 
we predict. Although preliminary optical characterizations 
have been performed"*^'"*^, linear dichroism has not yet been 
explored. 

Metliods 

DFT calculation. DFT calculations were performed using the generalized gradient 
approximation for the exchange-correlation potential, the projector augmented 
wave method^^'^^ and a plane-wave basis set as implemented in the Vienna ab- 
initio simulation package^^. The energy cutoff for the plane-wave basis was set to 
500 eV for all calculations. Two A:-meshes of 10 x 8 x 4 and 10 x 8 x 1 were 
adopted to sample the first BZ of the conventional unit cell of bulk and few-layer 
BP, and the mesh density of k points was kept fixed when calculating 
bandstructures using primitive cells. In optimizing the system geometry, van der 
Waals interactions were considered by the vdW-DF level with the optB88 exchange 
functional (optB88-vdW)^^'^'^. The shape and volume of each supercell were 
optimized fully and all atoms in the supercell were allowed to relax until the 
residual force per atom was < 0.001 eVA~^ Electronic bandstructures were 
calculated by the mBJ^^'^^ and hybrid functional (HSE06)^^'^^ methods based on 
the atomic structures obtained from the full optimization by optB88-vdW. The 
charge gradient used in the mBJ calculations, c= 1.1574, was extracted from the 
bulk value. 

Carrier mobility calculation. In 2D the carrier mobility is given by the expres- 
sion40-42 



where m* is the effective mass in the transport direction and m^is the average 
effective mass determined by ma = y^m*m*. The term Ei represents the 
deformation potential constant of the valence-band minimum for hole or 
conduction-band maximum for electron along the transport direction, defined by 
E\ = AVi/{Al/lo). Here AV, is the energy change of the i^^ band under proper cell 
compression and dilatation (calculated using a step of 0.5%), Iq is the lattice 
constant in the transport direction and A/ is the deformation of Iq. The elastic 
modulus C2D of the longitudinal strain in the propagation directions (both x and y) 
of the longitudinal acoustic wave is derived from {E — Eq)/So = C(A///o)^/2, where E 
is the total energy and So is the lattice volume at equilibrium for a 2D system. All 
structural properties in the calculation of carrier mobilities were obtained from 
optB88-vdW and properties related to the electronic structure were computed with 
the HSE06 functional. The temperature used for the mobility calculations was 
300 K. 

Uncertainty of deformation potential. In calculations of the electron and hole 
mobilities, the deformation potential is derived from linear fits to the respective 
energies of the conduction-band maximum and the valence-band minimum as 
functions of the lattice dilation or compression. Typical fitting results for E^y 
(electron state) and Ei^ (hole state) are shown in Supplementary Fig. 5a-c. 
The slope represents the deformation potential, which usually has an unavoidable 
uncertainty accounted for in the standard fitting error. Supplementary Figure 5d 
summarizes the relative errors for the deformation potentials of electrons and 
holes in few-layer BP for the x and y directions. For a given direction and type of 
carrier, the error is generally smallest for monolayer BP, while that for five-layer BP 
is the largest, except for the hole deformation potential along y. The error in 
Ely for hole states appears rather larger than the others, especially for monolayer 
BP where it exceeds 20%, presumably due to the small absolute value of Eiy 
for this case. Supplementary Fig. 5c shows details of the fit. The absolute value of 
Ely for hole states is only 0.15 eV, an order of magnitude smaller than the 
other energies; an error of 1 meV, believed to be the accuracy limit of DFT, 
may then lead to a fitting error of more than 10%. In this respect, a fitting 
error of 20% is definitely reasonable, and suggests that our calculation is highly 
accurate. 

Absorption spectra calculation. Absorption spectra were calculated from the 
dielectric function using the expression A{(jo) = 1 — e~ '^^'^^ ' where a{<jo)=^ is 

the absorption coefficient, n = y^ V^i + ^ ^^le index of refraction, Si and £2 are 
the real and imaginary parts of the dielectric function, co is the light frequency, c is 
the speed of light in vacuo and Az represents the unit-cell size in the c direction. 
The electronic structures were obtained from HSE06 results and the /c-mesh was 
doubled in calculating dielectric functions. Excitonic contributions were not con- 
sidered in our calculations. The total number of bands considered was set to be 
twice that used in the total-energy and bandstructure calculations. Because the 
dielectric function is a tensor, the absorption spectra along the three directions a 
(x), b iy) and c (z) were obtained separately. 
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